library(ggplot2)
library(pander)
library(WaveletComp)
setwd("~/Google Drive/Plodia_populations/Analysis") #NB set your working directory to whatever is appropriate
adults<-read.csv("Dead adult counts.csv")
adults$temp<-factor(adults$temp)
larvae<-read.csv("Larval counts.csv")
larvae$temp<-factor(larvae$temp)
larvae$food<-relevel(larvae$food, ref="good")
larvae$X5th<-larvae$X5F+larvae$X5M ## Combine male and female 5ths to give total
larvae$Total<-larvae$X3rd+larvae$X4th+larvae$X5th ## A few Totals don't equal the sum of 3rds, 4ths, 5ths
larvae2<-subset(larvae, number>10)
NB this analysis and the remaining analyses of population data are carried out using truncated time series with the first 10 weeks removed so as to avoid bias from transients associated with the initial setup of the populations.
Figure S1 shows time series for each population. All populations at 27 and 30 degrees persisted for the whole 82 weeks, as did two of the 33 degree populations, both on poor food. The populations that went extinct did so in less than 6 months, so the high temperature populations either went extinct quite rapidly or persisted for a long time.
rep<-ifelse(adults$box<7,1,ifelse(adults$box>12,3,2))
adults$food<-relevel(adults$food, ref="good")
adults<-data.frame(adults,rep)
rm(rep)
adults2<-subset(adults, week>10 & week<82)
temp_labels <-c("27" = "27°", "30" = "30°", "33" = "33°")
food_labels <-c("good" = "Standard diet", "bad" = "Poor diet")
rep_labels <-c("1" = "Replicate 1", "2" = "Replicate 2", "3" = "Replicate 3")
popplot<-ggplot(data=adults2, aes(x=week, y=total))
popplot<-popplot + geom_line(size = 0.2) + xlab("Time (weeks)") + ylab ("Total adults")
popplot<-popplot+facet_grid(rep~temp+food,
labeller = labeller(rep = rep_labels, temp = temp_labels, food =
food_labels))+theme_bw()+theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
strip.background =element_rect(fill="white", colour = "grey80"))
popplot
Figure S1 Time series for adult counts for each population. The first 10 weeks of data are not shown.
Figure S2 shows the mean population sizes (as weekly counts of dead adults) for each population. There are very clear differences between the temperatures, with the 30 degree treatments having the highest numbers of adult moths. The effect of food quality appears to be much smaller, if there is an effect at all. Of the 33 degree treatment populations, the two poor-quality food populations which persisted for the duration of the experiment can be seen to have much lower population sizes than the populations maintained at lower temperatures. The one good-quality population which has a mean population size which is similar to the two which persisted is a population which lasted longer than the other high-temperature ones but which eventually went extinct at around 26 weeks.
attach(adults2)
Mean.sizes<-tapply(total,box,function(x) mean(x, na.rm=TRUE))
Temperature<-factor(temp[1:18])
Food<-factor(food[1:18])
Sizes<-data.frame(Temperature, Food, Mean.sizes)
persisted<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)
Sizes2<-Sizes[persisted,]
sizeplot <- ggplot(Sizes2, aes(x=Temperature, y=Mean.sizes, fill=Food)) + ylim(0, 165)
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center",
position=position_dodge(0.5))
sizeplot <- sizeplot + theme_bw() + scale_fill_brewer(palette="Dark2",
name = "Food\nquality",
breaks = c("good", "bad"),
labels = c("Standard", "Poor"))
sizeplot<-sizeplot + ylab("Mean number of adults per week") + xlab ("Temperature °C")
sizeplot
rm(Mean.sizes, Temperature, Food)
detach(adults2)
Figure S2 Mean population sizes for the Plodia populations by temperature and food quality. The first 10 data points for each population have been removed to avoid transients.
We can fit a model to these data to assess effect sizes and the statistical significance of the effects of temperature and food quality. Because most of the high-temperature populations did not persist we’ll restrict the analysis to the 27 and 30 degree treatments only.
Full model
model1<-lm(Mean.sizes~Temperature*Food,
data=Sizes, subset=Temperature=="27" | Temperature=="30")
anova(model1)
Analysis of Variance Table
Response: Mean.sizes
Df Sum Sq Mean Sq F value Pr(>F)
Temperature 1 7851.2 7851.2 42.5928 0.000183 ***
Food 1 6.3 6.3 0.0344 0.857552
Temperature:Food 1 128.4 128.4 0.6966 0.428150
Residuals 8 1474.7 184.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(model1)
Minimal model
Model selection by removal of non-signficant terms leaves us with only temperature as an explanatory variable
model2<-lm(Mean.sizes~Temperature,
data=Sizes, subset=Temperature=="27" | Temperature=="30")
anova(model2)
Analysis of Variance Table
Response: Mean.sizes
Df Sum Sq Mean Sq F value Pr(>F)
Temperature 1 7851.2 7851.2 48.784 3.786e-05 ***
Residuals 10 1609.4 160.9
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rm(model2)
Since food quality has no significant effect we have a few more df to play with so we can include the two persistent populations at 33 degrees as well.
extinct<-c(5,11,17,18)
model3<-lm(Mean.sizes[-extinct]~Temperature[-extinct], data=Sizes)
summary(model3)
Call:
lm(formula = Mean.sizes[-extinct] ~ Temperature[-extinct], data = Sizes)
Residuals:
Min 1Q Median 3Q Max
-29.127 -6.317 3.265 6.449 13.183
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 93.322 4.963 18.802 1.04e-09 ***
Temperature[-extinct]30 51.157 7.019 7.288 1.57e-05 ***
Temperature[-extinct]33 -64.751 9.927 -6.523 4.29e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.16 on 11 degrees of freedom
Multiple R-squared: 0.9307, Adjusted R-squared: 0.9181
F-statistic: 73.87 on 2 and 11 DF, p-value: 4.206e-07
rm(model3)
Overall, temperature has a substantial effect on the mean population sizes. Increasing the temperature from 27 to 30 degrees increases mean population size roughly 1.5 times, but going from 27 degrees to 33 degrees means that those populations manage to persist at the higher temperature have a mean size which is roughly a third of the 27 degree populations and a fifth of the 30 degree ones.
attach(adults2)
sds<-tapply(total,box,function(x) sd(x, na.rm=TRUE))
Covars<-data.frame(Sizes, sds)
Covars$CV<-Covars$sds/Covars$Mean.sizes
persisted<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)
Covars<-Covars[persisted,]
sizeplot <- ggplot(Covars, aes(x=Temperature, y=CV, fill=Food))
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center",
position=position_dodge(0.5)) + theme_bw() + scale_fill_brewer(palette="Dark2",
name = "Food\nquality", breaks = c("good", "bad"),
labels = c("Standard", "Poor"))
sizeplot<-sizeplot + xlab("Temperature (°C)")
sizeplot <- sizeplot + ylab("Coefficient of variation for each adult population")
sizeplot
detach(adults2)
Figure S3 Coefficients of variation for each population
30º Populations were slightly less variable than 27º populations. The 33º populations that persisted had much higher CVs but this is possibly more of an effect of their low mean population size than anything else.
Can we discern any statistically significant patterns in the variance of the adult populations?
model1<-lm(CV~Temperature*Food, data=Covars, subset=Temperature=="27" | Temperature=="30")
drop1(model1, test="F")
Single term deletions
Model:
CV ~ Temperature * Food
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.0099351 -77.159
Temperature:Food 1 0.0014022 0.0113373 -77.575 1.1291 0.319
model2<-lm(CV~Temperature+Food, data=Covars, subset=Temperature=="27" | Temperature=="30")
drop1(model2, test="F")
Single term deletions
Model:
CV ~ Temperature + Food
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.011337 -77.575
Temperature 1 0.0090049 0.020342 -72.560 7.1485 0.02547 *
Food 1 0.0011779 0.012515 -78.389 0.9350 0.35882
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
No signficant interaction term and no effect of food. We can include the 33º populations in a final model.
model3<-lm(CV~Temperature, data=Covars)
drop1(model3, test="F")
Single term deletions
Model:
CV ~ Temperature
Df Sum of Sq RSS AIC F value Pr(>F)
<none> 0.013291 -91.436
Temperature 2 0.19094 0.204228 -57.186 79.014 2.978e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model3)
Call:
lm(formula = CV ~ Temperature, data = Covars)
Residuals:
Min 1Q Median 3Q Max
-0.049353 -0.021899 -0.007172 0.018844 0.056647
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.53089 0.01419 37.41 6.01e-13 ***
Temperature30 -0.05479 0.02007 -2.73 0.0196 *
Temperature33 0.29838 0.02838 10.51 4.47e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.03476 on 11 degrees of freedom
Multiple R-squared: 0.9349, Adjusted R-squared: 0.9231
F-statistic: 79.01 on 2 and 11 DF, p-value: 2.978e-07
CVs are significantly lower in the 30º populations and then significantly raised in the 33º populations. There is no discernable effect of food quality.
Figure S4 shows the mean counts per week for the total number of larvae per food section. By contrast with the data for adults there were fewer larvae in the 30 degree treatment than in the 27 degree one. The two 33 degree treatments which persisted havd very low numbers.
attach(larvae2)
Mean.sizes<-tapply(Total,box,function(x) mean(x, na.rm=TRUE))
Temperature<-factor(temp[1:18])
Food<-factor(food[1:18])
detach(larvae2)
Sizes<-data.frame(Temperature, Food, Mean.sizes)
Sizes<-Sizes[persisted,]
sizeplot <- ggplot(Sizes, aes(x=Temperature, y=Mean.sizes, fill=Food))
sizeplot <- sizeplot + geom_dotplot(binaxis="y", stackdir="center",
position=position_dodge(0.5)) + theme_bw() + scale_fill_brewer(palette="Dark2",
name = "Food\nquality", breaks = c("good", "bad"),
labels = c("Standard", "Poor"))
sizeplot <- sizeplot + ylab("Mean number of larvae per week")
sizeplot <- sizeplot + xlab("Temperature °C") + ylim(c(0,17))
sizeplot
Figure S4 Mean number of larvae counted per week per population - compare with Figure S2
Analysis of these data is similar to that for the adult populations. An initial model fitted to the 27 and 30 degree data only is reduced to one with only temperature as an explanatory variable and we can add in the two persisting 33 degree populations to that to give a final minimal adequate model with temperature as an explanatory variable but not food quality. By contrast with the adult populations, larval populations are smaller in the 30 degree treatment than in the 27 degree one, but as with the adult populations the two 33 degree populations which persisted have very low larval counts.
mod1<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)]*Food[-c(6,12)], data=Sizes)
anova(mod1)
mod2<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)]+Food[-c(6,12)], data=Sizes)
anova(mod2)
mod3<-lm(Mean.sizes[-c(6,12)]~Temperature[-c(6,12)], data=Sizes)
anova(mod3)
Minimal model
mod4<-lm(Mean.sizes~Temperature, data=Sizes)
anova(mod4)
Analysis of Variance Table
Response: Mean.sizes
Df Sum Sq Mean Sq F value Pr(>F)
Temperature 2 182.679 91.340 48.464 3.511e-06 ***
Residuals 11 20.732 1.885
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod4)
Call:
lm(formula = Mean.sizes ~ Temperature, data = Sizes)
Residuals:
Min 1Q Median 3Q Max
-1.73718 -0.79327 -0.05288 0.46554 3.14744
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.2564 0.5605 23.653 8.78e-11 ***
Temperature30 -2.6955 0.7926 -3.401 0.00592 **
Temperature33 -11.0353 1.1209 -9.845 8.64e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.373 on 11 degrees of freedom
Multiple R-squared: 0.8981, Adjusted R-squared: 0.8795
F-statistic: 48.46 on 2 and 11 DF, p-value: 3.511e-06
rm(Mean.sizes, Temperature, Food, persisted, mod1, mod2, mod3, mod4, Sizes)
Time averaged proportions of 3rd, 4th and 5th instar larvae are shown in Figure S4. There are clear differences between treatments, with the highest proportions of 5th instar larvae in the 27 degree good food populations and the lowest proportions being in the 30 degree populations.
prop.3rd<-larvae2$X3rd/larvae2$Total
prop.4th<-larvae2$X4th/larvae2$Total
prop.5th<-larvae2$X5th/larvae2$Total
#Calculate mean proportions per box
props.3rds<-tapply(prop.3rd,larvae2$box,function(x) mean(x, na.rm=TRUE))
props.4ths<-tapply(prop.4th,larvae2$box,function(x) mean(x, na.rm=TRUE))
props.5ths<-tapply(prop.5th,larvae2$box,function(x) mean(x, na.rm=TRUE))
ci.3rds.prop <- tapply(prop.3rd,
larvae2$box, function(x) qt(0.975,
df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))
#calculate 95%CI for each mean proportion per box
ci.4ths.prop <- tapply(prop.4th,
larvae2$box, function(x) qt(0.975,
df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))
ci.5ths.prop <- tapply(prop.5th,
larvae2$box, function(x) qt(0.975,
df = length(x) -1) * sd(x, na.rm = TRUE) / sqrt(length(x)))
ordered<-c(1,7,13,2,8,14,3,9,15,4,10,16,6,12)
props<-c(props.3rds[ordered],props.4ths[ordered],props.5ths[ordered])
cis <- c(ci.3rds.prop[ordered], ci.4ths.prop[ordered], ci.5ths.prop[ordered])
Instar<-factor(rep(c("3rd","4th","5th"),each=14))
Population<-factor(rep(c(1,2,3,1,2,3,1,2,3,1,2,3,1,2), times=3))
Food<-factor(rep(c("Standard diet","Standard diet","Standard diet",
"Poor diet","Poor diet","Poor diet","Standard diet","Standard diet",
"Standard diet","Poor diet","Poor diet","Poor diet","Poor diet",
"Poor diet"), times=3))
Food<-relevel(Food, ref = "Standard diet")
Temperature<-factor(rep(c(rep(c("27°","30°"), each=6),"33°","33°"), times=3))
Proportions<-data.frame(props,cis,Instar,Population,Food,Temperature)
## Comment out this line for all instars from 3rd to 5th
Proportions <- Proportions[Instar != "4th",]
p1<- ggplot(data=Proportions, aes(x=Population, y=props, fill=Instar))
## Remove position = "dodge" for a stacked plot
p1<- p1+geom_bar(stat="identity", position = "dodge") +ylim(0, 0.6)
#p1 <- p1 + geom_errorbar(aes(ymin = props - cis,
#ymax = props + cis), width = 0.1, position = position_dodge(0.9))
#Comment this out for a stacked plot
p1<- p1+facet_grid(~Temperature+Food,scales = "free_x",
space="free_x")+theme_bw() + ylab("Proportion of total larvae")
p1 <- p1 + xlab("Replicate") + theme(strip.background=element_rect( fill="white"))
#p1<- p1+scale_y_continuous(expand = c(0,0))+
#p1<- p1 + scale_fill_discrete(guide = guide_legend(reverse=TRUE))+scale_fill_brewer()
# test_palette <- c("chocolate3", "steelblue")
test_palette <- c("#ef8a62", "#67a9cf")
# p1<- p1 + scale_fill_brewer(guide = guide_legend(reverse=TRUE))
p1<- p1 + scale_fill_manual(values = test_palette)
p1
Figure S5 Proportions of the total larvae counted in each population that were in the 3rd or 5th instar. 4th instar proportions not shown for clarity.
These differences in the proportions are statistically significant. They can be analysed using a linear model with the mean number in one stadium (either 5th or 3rd instar larvae, it makes little difference) as a response variable and with the total number of larvae as a covariate.
mean.3rds<-tapply(larvae2$X3rd,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.4ths<-tapply(larvae2$X4th,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.5ths<-tapply(larvae2$X5th,larvae2$box,function(x) mean(x, na.rm=TRUE))
mean.all<-tapply(larvae2$Total,larvae2$box,function(x) mean(x, na.rm=TRUE))
keep<-c(1,2,3,4,7,8,9,10,13,14,15,16)
Vths<-mean.5ths[keep]
#Making variables with short names to keep the output readable
Total<-mean.all[keep]
#Subscripting with "keep" means that only the 27 and 30 degree boxes
#are being included - because there are only 2 33 degree ones and
#they are both poor food we can't test for an interaction if these are included.
Food<-larvae2$food[keep]
Temp<-as.factor(larvae2$temp[keep])
mod1<-lm(Vths~Total+Food*Temp)
anova(mod1)
Analysis of Variance Table
Response: Vths
Df Sum Sq Mean Sq F value Pr(>F)
Total 1 17.9399 17.9399 139.4633 7.079e-06 ***
Food 1 1.0627 1.0627 8.2611 0.02385 *
Temp 1 1.0342 1.0342 8.0401 0.02521 *
Food:Temp 1 0.8885 0.8885 6.9069 0.03401 *
Residuals 7 0.9004 0.1286
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod1)
Call:
lm(formula = Vths ~ Total + Food * Temp)
Residuals:
Min 1Q Median 3Q Max
-0.3994 -0.2969 0.0303 0.2409 0.3688
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9393 1.1555 0.813 0.44302
Total 0.4384 0.0831 5.275 0.00115 **
Foodbad -1.2829 0.3012 -4.260 0.00375 **
Temp30 -1.4403 0.3754 -3.836 0.00640 **
Foodbad:Temp30 1.0899 0.4147 2.628 0.03401 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3587 on 7 degrees of freedom
Multiple R-squared: 0.9587, Adjusted R-squared: 0.9352
F-statistic: 40.67 on 4 and 7 DF, p-value: 6.212e-05
rm(Vths,Total,Food,Temp)
Alternatively, an analysis using the proportion of 5th instar larvae as the response:
propmod<-lm(props~Food*Temperature, data=Proportions, subset=Instar=="5th")
anova(propmod)
Analysis of Variance Table
Response: props
Df Sum Sq Mean Sq F value Pr(>F)
Food 1 0.0153438 0.0153438 10.0619 0.011329 *
Temperature 2 0.0270404 0.0135202 8.8661 0.007455 **
Food:Temperature 1 0.0010241 0.0010241 0.6716 0.433659
Residuals 9 0.0137244 0.0015249
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
propmod2<-update(propmod,~.-Food:Temperature)
anova(propmod2)
Analysis of Variance Table
Response: props
Df Sum Sq Mean Sq F value Pr(>F)
Food 1 0.015344 0.0153438 10.4036 0.009091 **
Temperature 2 0.027040 0.0135202 9.1671 0.005476 **
Residuals 10 0.014749 0.0014749
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(propmod2)
Call:
lm(formula = props ~ Food + Temperature, data = Proportions,
subset = Instar == "5th")
Residuals:
Min 1Q Median 3Q Max
-0.047401 -0.027360 -0.002382 0.028821 0.046099
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.514059 0.019202 26.771 1.22e-10 ***
FoodPoor diet -0.079812 0.022172 -3.600 0.00485 **
Temperature30° -0.087631 0.022172 -3.952 0.00272 **
Temperature33° 0.007841 0.033259 0.236 0.81837
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0384 on 10 degrees of freedom
Multiple R-squared: 0.7419, Adjusted R-squared: 0.6644
F-statistic: 9.579 on 3 and 10 DF, p-value: 0.002749
The age structures of the larval populations are therefore significantly different between treatments, with a significant interaction between temperature and food quality. The 27 degree, good food populations have the most 5th instar larvae, and either increasing the temperature or changing to poor food decreases the relative abundance of 5th instar larvae in favour of younger animals. Food quality has little effect on age structure for the 30 degree populations, however, all of which have roughly similar age structures.
#### Plotting spectral densities with bootstrap 95% limits derived from the 95th
#quantiles for spectral densities of 10,000 replicates of white noise drawn from
#a negative binomial distribution with the mean and shape parameter equal to
#those for the time series in question
library(MASS)
adults<-read.csv("Dead adult counts.csv")
adults$temp<-factor(adults$temp)
rep<-ifelse(adults$box<7,1,ifelse(adults$box>12,3,2))
adults$food<-relevel(adults$food, ref="good")
adults<-data.frame(adults,rep)
rm(rep)
adults2<-subset(adults, week>10 & week<82)
attach(adults2)
sds<-tapply(total,box,function(x) sd(x, na.rm=TRUE))
means<-tapply(total,box,function(x) mean(x, na.rm=TRUE))
par(mfrow = c(3,4))
par(mar = c(1,1,1,1))
### Generate output data frame
Temperature <- c(rep(c(rep(c(27,30), each = 72),
rep(33, times = 36)), times = 2),
rep(c(27,30), each = 72))
Food <- c(rep(rep(c("Standard", "Poor", "Standard", "Poor", "Poor"),
each = 36), times = 2),
rep(c("Standard", "Poor", "Standard", "Poor"), each = 36))
Replicate <- c(rep(1, times = 180), rep(2, times = 180), rep(3, times = 144))
Freq <- numeric(0)
Spec <- numeric(0)
Quant <- numeric(0)
boxes<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)
for(i in boxes) {
box_test<-i
shape<-theta.md(adults2$total[adults2$box==box_test], mu = means[box_test],
dfr = length(adults2$total[adults2$box==box_test]-1))
# temp2<-rnorm(71000, mean = means[2], sd = sds[2])
temp2<-rnbinom(710000, mu = means[box_test], size = shape)
sims<-matrix(data = ifelse(temp2>=0, temp2, 0), nrow = 10000, ncol = 71)
spec1<-function(x, smoother = 3) {
smooth.spec = smoother
temp<-spectrum(sqrt(abs(x)), span = smooth.spec, plot = FALSE)
return(temp$spec)
}
spectra <- apply(sims, 1, spec1)
quant1 <- apply(spectra, 1, function(x) quantile(x, probs = c(0.95)))
smooth.spec<-3
temp1<-spectrum(sqrt(total[box==box_test & week <82]), span=smooth.spec, plot = FALSE)
Freq <- c(Freq, temp1$freq)
Spec <- c(Spec, temp1$spec)
Quant <- c(Quant, quant1)
# plot(temp1$freq, temp1$spec, type = "l")
#
# points(temp1$freq, quant1, type ="l", lty = 2)
#
# abline(v=1/6, lty=3, col="steelblue")
}
detach(adults2)
Spectral_output<-data.frame(Temperature, Food, Replicate, Freq, Spec, Quant)
Spectral_output$Temperature <- factor(Spectral_output$Temperature)
Spectral_output$Food <- factor(Spectral_output$Food)
Spectral_output$Food <- relevel(Spectral_output$Food, ref = "Standard")
Spectral_output$replicate <- factor(Spectral_output$Replicate)
par(mfrow = c(1,1))
temp_labels <-c("27" = "27°", "30" = "30°", "33" = "33°")
food_labels <-c("Standard" = "Standard food", "Poor" = "Poor food")
rep_labels <-c("1" = "Replicate 1", "2" = "Replicate 2", "3" = "Replicate 3")
p1 <- ggplot(data = Spectral_output, aes(x = Freq,
y = Spec)) + geom_line(size = 0.4) + theme_bw()
p1 <- p1 + geom_vline(xintercept = 1/6, colour = "steelblue",
linetype = "dotted", size = 0.5)
p1 <- p1 + geom_line(aes(x = Freq, y = Quant), linetype = "dashed", size = 0.4)
p1 <- p1 + xlab ("Frequency") + ylab("Spectral density")
p1 <- p1 + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),
strip.background =element_rect(fill="white", colour = "grey80"))
p1<-p1+facet_grid(Replicate ~ Temperature + Food,
labeller = labeller(Replicate = rep_labels,
Temperature = temp_labels, Food = food_labels),
scales = "free_x", space="free_x")
p1
Figure S6 Spectral densities for each population cage. Only those which persisted for the full period are shown. Dashed lines indicate approximate statistical significance levels: these show the 95th quantile for spectral densities of 10,000 replicates of white noise drawn from a negative binomial distribution with the mean and shape parameter equal to those for the time series in question. Dotted vertical lines indicate the frequency equivalent to a wavelength of six weeks. All populations at 27ºC and 30°C show strong peaks at a wavelength of six weeks, indicating regular cycles of one generation in length, apart from the 30°C populations with standard diet, which all show multiple, marginally significant peaks with little consistency. One 33°C population shows generation cycles but the other does not.
The data were square root transformed before analysis in order to normalise the error distributions. These smoothed spectral density plots show that in fact there are two important periodicities in the data - firstly a long-period cycle of between 35 and 25 weeks, as indicated by the peak on the left hand side of the graphs, and secondly short-period generation cycles of six weeks overlaid onto the longer term cycles. This is shown by the peak that in most cases aligns with the vertical dotted line, which indicates a frequency of 1/6 on each graph.
The locations of the peaks corresponding to the long-period cycles are given in table 2. Some populations show longer period cycles (approx. 35 weeks) and some shorter (approx 25 weeks). Note that because of the length of these cycles the spectral analysis will only fit peaks corresponding to periods of 72,36,24,18 etc. so the precise lengths of the long period cycles are hard to estimate.
popplot<-ggplot(data=adults2, aes(x=week, y=sqrt(total))) + labs(x="Time (weeks)",
y="Square root number of dead adults" )
popplot<-popplot + geom_line(size = 0.25)
popplot<-popplot+facet_grid(rep~temp+food)+geom_smooth(method="loess",
span=0.4, size=0.25, se=FALSE)
popplot<-popplot+theme_bw()+theme(panel.grid.major=element_blank(),
panel.grid.minor=element_blank())
popplot
Figure S7 The 18 population time series, square root transformed, with a fitted Loess smoother which makes the long-period cycles clear.
There is no apparent association between temperature and cycle length, but there is between food quality and cycle length. If we consider the 27 and 30 degree populations, all the poor food populations show short cycles (six out of six), whereas most of the good food populations show longer cycles (4 out of 6). This is not quite statistically significant:
cycles<-matrix(c(6,0,2,4), ncol=2)
rownames(cycles)<-c("Long.cycles","Short.cycles")
colnames(cycles)<-c("Poor.food","Good.food")
cycles
Poor.food Good.food
Long.cycles 6 2
Short.cycles 0 4
fisher.test(cycles)
Fisher's Exact Test for Count Data
data: cycles
p-value = 0.06061
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
0.9097849 Inf
sample estimates:
odds ratio
Inf
Regarding the generation cycles, the spectral densities reveal strong peaks corresponding to a period of 6 weeks in all of the 27 degree populations - these populations have dynamics that are essentially dominated by a long period cycle with the shorter period generation cycles overlaid onto them, and there is little indication of other important periodicity. For the 30 degree populations, there are strong peaks at six weeks in all three of the “poor food” populations although the picture is not as simple for these populations and population 14 has a notable third peak corresponding to a period of roughly 10 weeks. The 30 degree populations given “good” food have complex spectra with multiple peaks and little indication of strong generation cycles.
Of the 33 degree populations which persisted, one (population 6) had little to indicate generation-length cycles, whereas the other (population 12) had a strong peak at a frequency corresponding to a six-week period.
Wavelet analysis (Cazelles et al. 2008) allows us to examine whether any periodic behaviour exhibited by our populations changed during the 82 weeks they were maintained, and is used here to try to understand the differences in periodic behaviour between the 27ºC and 30ºC populations. Because the long-period cycles in these populations have wavelengths which are greater than 25% of the total length of the time series, these were removed from the data by smoothing with a loess function of span 0.4 prior to analysis. Wavelet transformations were carried out using the WaveletComp package (Roesch & Schmidbauer 2014) using a Morlet wavelet and plotted with 5% significance contours calculated using the default (white noise) options.
library(WaveletComp)
######## Plotting the wavelets without the population data.
par(mfrow = c(3,5))
par(mar=c(1,1,2,1)) #Margin sizes
to.plot<-c(1,2,3,4,6,7,8,9,10,12,13,14,15,16)
#Only plot the 27 and 30 degree populations plus 2 33 degree ones which persisted
for (i in to.plot) {
box<-adults2$total[adults2$box==i]
weeks<-adults2$week[adults2$box==i]
temp<-loess(box~weeks, span=0.4)
#plot(temp$residuals, type="l", main=paste("Population ", i), xaxs="i")
#Plot the residuals after a loess fit -
#in other words the data with the long wavelength cycles removed
mydata<-data.frame(x=box)
my.w2 = analyze.wavelet(mydata, "x",
#Calculate wavelet power spectrum
loess.span = 0.4,
#Remove long wavelength cycles before analysis
dt = 1, dj = 1/250,
lowerPeriod = 2,
#Check for periods between 2 and 24 weeks
upperPeriod = 24,
make.pval = T, n.sim = 10)
#Use more simulations for the significance test
#than the default which gives smoother contours
wt.image(my.w2, color.key = "interval",
#Plot wavelets on an interval scale
n.levels = 250,color.palette = "rainbow(n.levels, start=0, end=0.7)",
legend.params = list(lab = "wavelet power levels", mar = 4.7),
graphics.reset=FALSE,
#Need this so that it doesn't just overplot in the same place
plot.ridge=FALSE,
#Plot "Ridge of power"
lvl=0.8,
#Only plot ridge for reasonably high powers
siglvl=0.01,
#Set p-value for significance contour to 0.05 rather than the default 0.1
label.time.axis = FALSE, label.period.axis = FALSE,
plot.legend=FALSE)
mtext(text = paste("Temp = ", adults2$temp[adults2$box == i][1],
" Food = ", adults2$food[adults2$box == i][1]), cex = 1.3)
#Add text above plot giveing temperature and food quality
}
par(mfrow=c(1,1))
Figure S8. Wavelet power spectra of population time series. The x-axis gives time in weeks and the y-axis indicates the period. The power is indicated by colour, with blue indicating low power whereas green, yellow and red show increasing power, and the paler regions at the left and right hand sides of the plots show the regions where effects associated with the beginning and end of the time series reduce certainty. The time series identified by spectral analysis (Fig. 4) as showing 6-week generation cycles, including all of the 27°C populations, the 30°C poor diet populations and the second replicate of the 33°C poor diet treatment all show areas of high power at a period of 6 weeks, with those with particularly strong 6-week peaks in the spectral analysis showing ‘ridges’ of power across the time series, indicating that generation cycling is occurring for most of the period that the populations were maintained. The 30°C populations maintained on standard diet, by contrast, show brief periods of regularity at a variety of periods, and only the second one shows any peak of power near the 6-week periodicity found in the other time series. The 33°C population which does not show generation cycles (the first one) has a very brief period of cycling at the start of the time series after which there is no real periodicity. Populations which became extinct during the experiment are not shown.